function [Cov0_dX_m_n_cond_on_bus_cycle,corr0_dX_m_n_cond_on_bus_cycle] = ...
    CalcTimeSeriesCorr1_cond_on_bus_cycle(IDX_RECESSION_THRESH,X_Nxz,prob_Nxz,Phi_Nxz,m,n,...
    Phi_Nxz_m,Phi_Nxz_n)
% computes 
% corr(X(t+m)-X(t),X(t+m+n)-X(t+m)|REC(t)) = A/(sqrt(B)*sqrt(C))
% 
% where REC(t) = z(t)<=z_thresh
% 
% A = cov(X(t+m)-X(t),X(t+m+n)-X(t+m)|z(t)=z)
%   = cov(X(t+m),X(t+m+n)|z(t)=z)-var(X(t+m)|z(t))
%     - cov(X(t),X(t+m+n)|z(t)) + cov(X(t),X(t+m)|z(t))
% B = var(X(t+m)-X(t)|z(t)=z)
%   = var(X(t+m)|z(t))+var(X(t)|z(t))-2*cov(X(t+m),X(t)|z(t))
% C = var(X(t+m+n)-X(t+m)|z(t)=z)
%   = var(X(t+m+n)|z(t))+var(X(t+m)|z(t))-2*cov(X(t+m+n),X(t+m)|z(t))

    [NN Nx Nz] = size(prob_Nxz);

    % recession and expansion ranges
    IDX_REC = 1:IDX_RECESSION_THRESH;
    IDX_EXP = IDX_RECESSION_THRESH+1:Nz;

    IDX_REGIME = cell(2,1);
    IDX_REGIME{1} = IDX_REC;
    IDX_REGIME{2} = IDX_EXP;

    if nargin<=6
        Phi_Nxz_m = mpower2(Phi_Nxz,m);
    end
    
    if nargin<=7
        Phi_Nxz_n = mpower2(Phi_Nxz,n);
    end
    
    Phi_Nxz_mn = Phi_Nxz_m*Phi_Nxz_n;
    
    prob_z = reshape(squeeze(sum(prob_Nxz,[1 2])),[],1);
    prob_regime = zeros(2,1);
    prob_regime(1) = sum(prob_z(IDX_REC));
    prob_regime(2) = sum(prob_z(IDX_EXP));
    
    % E[X(t)|REC(t)] and Var[X(t)|REC(t)]
    E0_X0 = zeros(2,1); Var0_X0 = zeros(2,1);
    for i = 1:2
        E0_X0(i) = sum(prob_Nxz(:,:,IDX_REGIME{i}).*X_Nxz(:,:,IDX_REGIME{i}),'all')/prob_regime(i);
        Var0_X0(i) = sum(prob_Nxz(:,:,IDX_REGIME{i}).*(X_Nxz(:,:,IDX_REGIME{i}).^2),'all')/prob_regime(i) ...
            - E0_X0(i)^2;
    end
    
    % E[X(t+m)|REC(t)] and Var[X(t+m)|REC(t)]
    E0_Xm = zeros(2,1); Var0_Xm = zeros(2,1);
    temp_m = reshape(Phi_Nxz_m*X_Nxz(:),[NN Nx Nz]); % E_t[X(t+m)]
    for i  = 1:2
        E0_Xm(i) = sum(prob_Nxz(:,:,IDX_REGIME{i}).*temp_m(:,:,IDX_REGIME{i}),'all')/prob_regime(i);
        Var0_Xm(i) = sum(prob_Nxz(:,:,IDX_REGIME{i}).*(temp_m(:,:,IDX_REGIME{i}).^2),'all')/prob_regime(i) ...
            - E0_Xm(i)^2;
    end
    
    % E[X(t+m+n)|REC(t)] and Var[X(t+m+n)|REC(t)]
    E0_Xmn = zeros(2,1); Var0_Xmn = zeros(2,1);
    temp_mn = reshape(Phi_Nxz_mn*X_Nxz(:),[NN Nx Nz]); % E_t[X(t+m)]
    for i  = 1:2
        E0_Xmn(i) = sum(prob_Nxz(:,:,IDX_REGIME{i}).*temp_mn(:,:,IDX_REGIME{i}),'all')/prob_regime(i);
        Var0_Xmn(i) = sum(prob_Nxz(:,:,IDX_REGIME{i}).*(temp_mn(:,:,IDX_REGIME{i}).^2),'all')/prob_regime(i) ...
            - E0_Xmn(i)^2;
    end
    
    % Cov(X(t),X(t+m)|REC(t))
    Cov0_X0_Xm = zeros(2,1);
    for i = 1:2
        Cov0_X0_Xm(i) = sum(prob_Nxz(:,:,IDX_REGIME{i}).*X_Nxz(:,:,IDX_REGIME{i}).*temp_m(:,:,IDX_REGIME{i}),'all')/prob_regime(i) ...
            - E0_X0(i)*E0_Xm(i);
    end
    
    % Cov(X(t),X(t+m+n)|z(t))
    Cov0_X0_Xmn = zeros(2,1);
    for i = 1:2
        Cov0_X0_Xmn(i) = sum(prob_Nxz(:,:,IDX_REGIME{i}).*X_Nxz(:,:,IDX_REGIME{i}).*temp_mn(:,:,IDX_REGIME{i}),'all')/prob_regime(i) ...
            - E0_X0(i)*E0_Xmn(i);
    end
    
    % Var[X(t+m)-X(t)|REC(t)]
    Var0_diff_Xm_X0 = Var0_Xm + Var0_X0 - 2*Cov0_X0_Xm;
    
    % Cov(X(t+m),X(t+m+n)|REC(t))
    Cov0_Xm_Xmn = zeros(2,1);
    temp_m_n = reshape(Phi_Nxz_m*(X_Nxz(:).*(Phi_Nxz_n*X_Nxz(:))),[NN Nx Nz]);
    for i = 1:2
        Cov0_Xm_Xmn(i) = sum(prob_Nxz(:,:,IDX_REGIME{i}).*temp_m_n(:,:,IDX_REGIME{i}),'all')/prob_regime(i) ...
            - E0_Xm(i)*E0_Xmn(i);
    end
    
    % Var[X(t+m+n)-X(t+m)|z(t)]
    Var0_diff_Xmn_Xm = Var0_Xmn + Var0_Xm- 2*Cov0_Xm_Xmn;
    
    Cov0_dX_m_n_cond_on_bus_cycle = Cov0_Xm_Xmn - Var0_Xm - Cov0_X0_Xmn + Cov0_X0_Xm;

    corr0_dX_m_n_cond_on_bus_cycle = Cov0_dX_m_n_cond_on_bus_cycle./(sqrt(Var0_diff_Xm_X0).*sqrt(Var0_diff_Xmn_Xm));

end